\chapter{Input parameters}\label{ch:inputParameters}
\index{file!input data}

This chapter describes the main input files for Brahms: {\tt brahms.md}, {\tt lipid.ff}, {\tt water.ff} and {\tt brahms.an}. Where applicable, units are reported between square brackets. 
Some parameters are marked as {\em logical}, meaning that they can be either TRUE (when set to 1 or any other positive integer) or FALSE (for any other values, including when left blank). They typically control the activation of various functionalities. The description given below refers to their behavior when set TRUE; when set FALSE, the corresponding functionality is simply switched off.  
There are currently no default values assigned by Brahms; if a parameter is not set, it is normally set to 0 by the compiler. 
Important things to note:
\begin{itemize}
\item the input files should contain no blank lines;
\item the last line of the input files should terminate with a ``newline'' character.
\end{itemize}
The examples provided in the \brahms package should clarify these rules. In general, it is obviously crucial that the input parameters are correctly loaded in the program. At the start of a simulation, \brahms loads all input parameters and subsequently writes them out to the file {\tt inputParameters.log}; in case the simulation behaves unexpectedly, {\tt inputParameters.log} should be checked against the input files to make sure the parameters have been properly loaded (see also section~\ref{sec:outFiles}).

Many of \brahms input parameters have been implemented following Rapaport;~\cite{rapa} therefore, to understand how they work, it is useful to look them up also in the original reference~\citep{rapa}[e.g., pp.~509-515]. For a really full (yet possibly onerous) understanding, parameters can be traced in the source code.\footnote{For example, try {\tt grep anyParameterName *.c *.h} on the command line from the \brahms source code ({\tt src/}) directory.} 

\section{{\tt brahms.md}}\label{sec:brahmsmd}
The input file {\tt brahms.md} contains general \underline{m}olecular \underline{d}ynamics parameters, together with parameters used to generate a system from scratch. See following paragraphs for individual descriptions.

\paragraph{{\tt deltaT}} [fs] Integration timestep. The standard value to use with \brahms is 15\,fs. Smaller time steps will yield better energy conservation but a lower sampling speed, possibly resulting in a waste of computer time. However, short time steps (even 0.1\,fs or shorter) can sometimes be required to equilibrate a (newly generated) system. Also, time steps of 1-10\,fs are sometimes necessary to run under ``extreme'' conditions (such as very high temperature). Timesteps larger than 15\,fs will improve sampling speed but might compromise stability and proper energy conservation.~\citep{winger09}

\paragraph{{\tt stepAvg}} [number of steps] Output frequency of text summary printed on screen. For development or checking purposes, any output frequency from 1 to 1000 steps is typically used. For long ``production'' jobs, 100000 steps or more are normally a better choice.
 
%\paragraph{{\tt stepEquil}} [number of steps] Duration of initial equilibration phase, characterized by temperature control and fixed volume. This is really only useful in simple example runs of pure water systems. Membrane simulation are typically equilibrated through standard runs in the NPT ensemble.

\paragraph{{\tt stepLimit}} [number of steps] Total simulation length. See section~\ref{sec:restarting} about possible ambiguities when restarting a simulation.

\paragraph{{\tt resetTime}} ({\em logical}) When restarting a simulation, resets the simulation time to 0.
 
\paragraph{{\tt doCheckpoint}} ({\em logical}) Activates checkpointing (Rapaport~\citep{rapa} pp.~500-504). This involves two features: i) if checkpoint files are available, \brahms will restart the corresponding simulation; ii) new checkpoint files will be generated and updated every {\tt stepCheckpoint} steps (see below).

\paragraph{{\tt stepCheckpoint}} [number of steps] Frequency at which checkpoint files are written (Rapaport~\citep{rapa} pp.~500-504).

\paragraph{{\tt stepPdb}} [number of steps] Frequency at which the system coordinates, in {\tt pdb} format, are written to the {\tt trajectory.pdb} output file. This file can be used to visualize the system (for instance using VMD~\citep{HUMP96}) or to carry out analysis. Warning: a low value of {\tt stepPdb} can lead to the generation of huge trajectory files. For example, for membrane simulations, since most of the analysis is carried out automatically  by \brahms while the program is running, a frequency of 100000 can be enough.

\paragraph{{\tt runId}} Simulation identifier (Rapaport~\citep{rapa} p.~501). 

\paragraph{{\tt applyThermostat}} ({\em logical}) Activates thermostatting; the temperature will be maintained at {\tt extTemperature} (see below). The temperature-control algorithm used in \brahms is Berendsen's weak-coupling scheme.~\citep{ber84} The statistical mechanical ensemble for the weak-coupling thermostat is intermediate between the canonical and the microcanonical ensembles.~\citep{morishita00} 

\paragraph{{\tt extTemperature}} [\textcelsius] ``External'' temperature of the thermostat. When the thermostat is on ({\tt applyThermostat} set to 1) the system temperature is maintained at this desired value.

\paragraph{{\tt tauT}} [fs] Temperature coupling constant of Berendsen's weak-coupling scheme.~\citep{ber84} Typically set in relation to {\tt deltaT}, that is, {\tt tauT} should normally be 10-30 times larger than {\tt deltaT}. For example, if {\tt deltaT} is set to 15\,fs, {\tt tauT} should be set to~$\approx200-500$\,fs. 

\paragraph{{\tt applyBarostat}} ({\em logical}) Activates barostatting; the pressure will be maintained at {\tt extPressure} (see below). The pressure-control algorithm used in \brahms is Berendsen's weak-coupling scheme.~\citep{ber84} 

\paragraph{{\tt extPressure}} [atm] ``External'' pressure of the barostat. When the barostat is on ({\tt applyBarostat} set to 1) the system pressure is maintained (on average) at this desired value. Note that it is normal for the {\em instantaneous} pressure to oscillate significantly, especially for small systems. The important thing is that the {\em average} value (over a sufficiently large number of time steps) is kept constant. 

\paragraph{{\tt tauP}} [fs] Pressure coupling constant of Berendsen's weak-coupling scheme.~\citep{ber84} Typically set in relation to {\tt deltaT}, that is, {\tt tauP} should normally be 20-60 times larger than {\tt deltaT}.

\paragraph{{\tt flexBox}} ({\em logical}) Allows box flexibility in response to pressure changes. This should normally be activated for bilayer simulations. Self-assembly runs might be more stable if {\tt flexBox} is set to 0.

\paragraph{{\tt keepTetragonal}} ({\em logical}) Keeps the box ``tetragonal'' - the box edges are constrained to remain perpendicular to each other. {\em This should always be activated when using the barostat}, because currently \brahms cannot deal with non-tetragonal boxes. This limitation is not a problem for most simulations of systems in the liquid phase, including fluid-phase bilayers.\footnote{However, non-tetragonal boxes can be more realistic when simulating solids.}

\paragraph{{\tt keepSquare}} ({\em logical}) Constrains the relative size of the x and y edges of the simulation box. Typically, this constrains the membrane interfacial plane (conventionally the xy plane) to remain a square. {\em This should normally be activated for pre-assembled membrane simulations}; it does not make them any less realistic (because of bilayers' ``natural'' homogeneity and symmetry in the xy plane), yet it prevents potential instabilities (caused by possible box ``thinning'' leading to violations of the minimum image convention).

\paragraph{{\tt rCutLipLip}} [nm] Cutoff radius for nonbonded interactions within lipid sites and between lipid and water sites. The recommended value for simulations of lipids modeled with the ELBA force field is 1.2.


\paragraph{{\tt rCutWatWat}} [nm] Cutoff radius for nonbonded interactions between water sites. The recommended value is 0.9.

\paragraph{{\tt rCutSolute}} [nm] Cutoff radius for nonbonded interactions between solute sites (not currently used).

\paragraph{{\tt rCutSoluteElse}} [nm] Cutoff radius for nonbonded interactions between solute sites and any other site (not currently used).

\paragraph{{\tt rNebrShell}} [nm] Thickness of the ``shell'' outside the cutoff sphere used to build the neighbor list (Rapaport~\citep{rapa} pp.~54-58). ``Typical'' values are 0.1-0.2; optimal values can be found for each application through trial-and-error refinement tests. This parameter can affect the accuracy and efficiency of the dynamics integration; see discussion below for the related parameter {\tt stepNebr}.

\paragraph{{\tt stepNebr}} [number of steps] Neighbor list update frequency. In other words, the neighbor list will be rebuilt every {\tt stepNebr} timesteps. Typical values (for a 15-fs timestep) are 5-20. In general, small values are more likely to guarantee that all interactions within cutoff range are properly considered, but the frequent reconstructions of the list will slow down the simulation. On the other hand, large values of {\tt stepNebr} will increase the simulation speed but might compromise stability and energy conservation. It is worth testing empirically which combination of {\tt rNebrShell} and {\tt stepNebr} yields the best compromise between accuracy (typically in terms of energy conservation) and efficiency. 

\paragraph{{\tt nebrTabFac}} Determines how much storage should be provided for the neighbor list, per site (Rapaport~\citep{rapa} pp.~54-58). Typically 100-300. This parameter should not affect simulation speed or indeed any result.

\paragraph{{\tt removeSystemTranslation}} ({\em logical}) Removes any net translation of the whole system at every step. This should be normally activated, to prevent potentially serious artefacts.~\citep{harve98a} However, this should be off when checking energy conservation in the NVE ensemble.

\paragraph{{\tt removeMonolayersTranslation}} ({\em logical}) In a simulation of a lipid bilayer, removes any net translation of each of the two monolayer (at every step). There is debate whether this operation is needed/justified.~\citep{patra04a,klauda06b,roark09} This removal should probably be activated for small systems ($<128$ lipids) and for gel-phase bilayers, while it should make no difference for large ($>288$ lipids) membranes in the fluid phase. For sure this should be off when checking energy conservation. Warning: this operation relies on the assumption that the bilayer plane is parallel to the xy plane defined by the coordinates' system; while this is normally the case for preassembled bilayer simulations in Brahms, self-assembled systems might not be oriented along the xy plane, in which case {\tt removeMonolayersTranslation} should definitely be off.

\paragraph{{\tt nSites}} To be set corresponding to the total number of sites in the simulated system.

\paragraph{{\tt nWaters}} To be set corresponding to the total number of water sites in the simulated system.

\paragraph{{\tt nDOPCsDSPCs}} To be set corresponding to the total number of DOPC or DSPC lipids in the simulated system. 

\paragraph{{\tt nDOPEs}} To be set corresponding to the total number of DOPE lipids in the simulated system. 

\paragraph{{\tt nLipids}} To be set corresponding to the total number of lipids in the simulated system. 

\paragraph{{\tt nSolutes}} To be set corresponding to the total number of solutes in the simulated system. 

\paragraph{{\tt nTypes}} To be set corresponding to the total number of site types in the simulated system. Currently there are 7 types in the ELBA force field.

\paragraph{{\tt region}} [nm nm nm] Used for the generation of a system from scratch - see section~\ref{sec:initLip}.

\paragraph{{\tt adjustRegion}}  ({\em logical}) Changes the box edges to {\tt regionAdjusted} (see below). Not normally used.

\paragraph{{\tt regionAdjusted}} [nm nm nm] Used to modify the simulation box edges when restarting a simulation. Not normally used.

\paragraph{{\tt randSeed}} Random number seed (Rapaport~\citep{rapa} pp.~491-492).

\paragraph{{\tt initHalfCellWat}}  Used for the generation of a membrane system from scratch - see section~\ref{sec:initLip}.

\paragraph{{\tt initUcell}} Used for the generation of a water system from scratch - see section~\ref{sec:initWat}.

\paragraph{{\tt loadStructure}} Not currently used.

\paragraph{{\tt loadVelocities}} Not currently used.

\paragraph{{\tt centerInputStruct}} Not currently used.

\paragraph{{\tt reCenterBilayer}}  ({\em logical}) ``Pulls'' gradually the bilayer to the origin of the z axis; this can be used after self-assembly simulations, which typically yield a bilayer which is not in the center of the simulation box.

\paragraph{{\tt zConstraint}} Not currently used.

\paragraph{{\tt insertSolute}} Not currently used.

\section{{\tt lipid.ff}}
The input file {\tt lipid.ff} contains the \underline{f}orce\underline{f}ield parameters of the ELBA model for coarse-grain lipids.~\cite{orsi11elba} This file is supplied with the \brahms package and should not be changed, unless the user really wants to change the force field.

\section{{\tt water.ff}}
The input file {\tt water.ff} contains the \underline{f}orce\underline{f}ield parameters of the ELBA model for coarse-grain water.~\cite{orsi11elba} This file is supplied with the \brahms package and should not be changed, unless the user really wants to change the force field.

\section{{\tt brahms.an}}\label{sec:brahms.an}
The optional input file {\tt brahms.an} controls various \underline{an}alysis operations that \brahms can carry out ``on the fly'', that is, while the simulation is running. All parameters of {\tt brahms.an} are described in the following; subsections collect related sets of parameters.

\subsection{Region dimensions}\label{sec:regDims} 

\paragraph{{\tt writeAreaVol}}  ({\em logical}) Write area and volume of the simulation region at every step on  file ({\tt area$\_$volume.dat}). At the end of the run, these data can be used to calculate the area and volume compressibility moduli.~\cite{marioDmpcDopc} A script to carry out this calculation is distributed with the \brahms package, in the directory {\tt analysisScripts/general}. Open the scripts with any text editor and follow the usage instructions at the top of the file.

\subsection{Lipid lateral diffusion} \label{sec:lld}
These parameters are related to the ``lateral'' lipid diffusion process, that is, the motion of the center of mass of each individual lipid molecule in the membrane plane (conventionally, the xy plane). To measure diffusion, \brahms implements the tools described by Rapaport~\cite{rapa}~(pp.~120-128); please check this reference for more details.
Note that the algorithms underlying the diffusion measurements in \brahms work only for uninterrupted simulations; restarted runs do not retain any memory from previous (interrupted) diffusion measurements.
 
\paragraph{{\tt latDiff}} ({\em logical}) Activates lipid lateral diffusion measurements.
\paragraph{{\tt nBuffLatDiff}} Number of sets of data being collected at any time.
\paragraph{{\tt nValLatDiff}} Number of measurements contributing to the set used to produce a single unaveraged estimate of the diffusion coefficient. 
\paragraph{{\tt stepLatDiff}} [number of steps] Measurement frequency.
\paragraph{{\tt limitLatDiffAv}} Total number of individual estimates used to produce an average value of the diffusion coefficient; such an average is written on file. 
\paragraph{{\tt writeLipLatMotion}} ({\em logical}) Writes on file the mass center coordinates of each lipid in the ``upper'' monolayer every {\tt stepLatDiff} steps. The file produced can be used to plot and analyze single lipid lateral traces projected onto the {\em xy} plane. For an example, see Fig.~6 in Orsi et al.~\cite{marioDmpcDopc} Scripts to produce such diagrams with the Xmgrace plotting program are distributed with the Brahms package - see {\tt analysisScripts/general/lipidDiffusionTraces.txt}. Note that it would not be possible to obtain the same results from the trajectory file, because this file contains coordinates that are ``wrapped'' around the periodic boundaries. Instead, {\tt writeLipLatMotion} generates and updates auxiliary structures which maintain ``unwrapped''  coordinates throughout the simulation. 

\subsection{Transmembrane profiles} 

This section describes the {\tt brahms.an} parameters controlling measurements of properties as a function of ``depth'' inside the membrane. The general procedure involves considering ``slices'' of the system along planes %parallel to the $xy-$plane, i. e. 
perpendicular to the $z$ axis (the interface normal by convention). 
Several bilayer properties are homogeneous inside a particular slice, due to the intrinsic axial symmetry of the system. %. Hence these properties depend only on the position along the normal, 
Therefore single curves, {\em profiles} evaluated as a function of $z$, provide full characterization. Typical membrane profiles are: electron density, lateral pressure, electric field, water polarization, and electrostatic potential.

The measurement process for all profiles involve similar parameters, here generically called {\tt sizeHistProfile}, {\tt stepProfile}, {\tt limitProfile}. Each profile is constructed as a histogram, with each histogram ``bin'' corresponding to a particular slice of the system. The parameter {\tt sizeHistProfile} sets the bin number (which corresponds to the number of ``slices''). More (fewer) bins yield an increased (decreased) measurement resolution. The relation between resolution, number of bins and system size is: $sizeHistProfile\times resolution = zEdge $, with $zEdge$ the length of the simulation box edge perpendicular to the membrane plane (by convention this is the z-dimension\footnote{For example, a desired resolution of 0.1\,nm for a box which measures 6.4\,nm in the $z$ dimension can be obtained by setting {\tt sizeHistProfile} to $zEdge / resolution = 6.4 / 0.1 = 64 $.}).   
The parameter {\tt stepProfile} sets the desired number of steps between single measurements.\footnote{For example, profiles can be evaluated at every step by setting {\tt stepProfile} to 1;  usually however this is not necessary, and in fact such a high frequency will slow down the calculations. More appropriate values of {\tt stepProfile} for ``long'' runs are 100-1000. The ``optimum'' value is to be found empirically through test runs.} The parameter {\tt limitProfile} sets the number of single measurements used by \brahms to calculate an average which is then written on file.\footnote{For example, take a 100\,ns production run, corresponding to 10\,000\,000 steps with a 10\,fs timestep. Setting {\tt stepProfile} to 10 and {\tt limitProfile} to 10\,000 will produce 100 output files, each representing an average over a 1-ns ``block'' (corresponding to $10\times10\,000=100\,000$\,steps). For the same run, setting {\tt stepProfile} to 100 and {\tt limitProfile} to 10\,000 will produce 10 output files, each representing an average over a 10-ns ``block'' (corresponding to $100\times10\,000=1\,000\,000$\,steps). Note that you would typically want the total average over the entire simulation; for a specific value of {\tt stepProfile}, such an average would not depend on {\tt limitProfile}.}    

\subsubsection{Electron density profiles} 
These parameters control the calculation of electron density profiles for the whole system and for individual site types; see Fig.~2 in Orsi et al.~\cite{marioDmpcDopc} and  Fig.~6 in Orsi et al.~\cite{mario08}

\paragraph{{\tt edp}} ({\em logical}) Activates profile calculation.
\paragraph{{\tt sizeHistEdp}} Resolution of the electron density profile histograms.
\paragraph{{\tt stepEdp}} [number of steps] Interval between individual profile evaluations (measurement frequency).
\paragraph{{\tt limitEdp}} Number of individual evaluations used to produce a single average curve which is then written on file.

\subsubsection{Electrostatic potential profiles}
These parameters control the calculation of electrostatic potential profiles for the whole system and for individual site types; see Fig.~4 in Orsi et al.~\cite{marioDmpcDopc} and  Fig.~10 in Orsi et al.~\cite{mario08}

\paragraph{{\tt epp}} ({\em logical}) Activates profile calculation.
\paragraph{{\tt sizeHistEpp}}  Resolution of the electron electrostatic potential histograms.
\paragraph{{\tt stepEpp}} [number of steps] Interval between individual profile evaluations (measurement frequency).
\paragraph{{\tt limitEpp}} Number of individual evaluations used to produce a single average curve which is then written on file.

\subsubsection{Lateral pressure profile}
These parameters control the calculation of the lateral pressure profile; see Fig.~3 in Orsi et al.~\cite{marioDmpcDopc} and  Fig.~7 in Orsi et al.~\cite{mario08}
\paragraph{{\tt lpp}} ({\em logical}) Activates profile calculation.
\paragraph{{\tt sizeHistLpp}} Resolution of the lateral pressure profile histograms.
\paragraph{{\tt stepLpp}} [number of steps] Interval between individual profile evaluations (measurement frequency).
\paragraph{{\tt limitLpp}} Number of individual evaluations used to produce a single average curve which is then written on file.

\subsubsection{Water polarization profile}
These parameters control the calculation of the water polarization profile; see Fig.~9 in Orsi et al.~\cite{mario08}
\paragraph{{\tt wpp}} ({\em logical}) Activates profile calculation.
\paragraph{{\tt sizeHistWpp}} Resolution of the water polarization profile histograms.
\paragraph{{\tt stepWpp}} [number of steps] Interval between individual profile evaluations (measurement frequency).
\paragraph{{\tt limitWpp}} Number of individual evaluations used to produce a single average curve which is then written on file.

\subsection{Radial distribution function in pure water systems}

\brahms implements the RDF measurement method from Rapaport~\cite{rapa}~(pp.~222-225).

\paragraph{{\tt rdfWat}} ({\em logical}) Activates RDF calculation.
\paragraph{{\tt rangeWatRdf}} [nm] Upper limit of RDF measurement.
\paragraph{{\tt sizeHistWatRdf}} Resolution of the RDF histograms.
\paragraph{{\tt stepWatRdf}} [number of steps] Interval between individual RDF evaluations.
\paragraph{{\tt limitWatRdf}} Number of individual evaluations used to produce a single average which is then written on file.

\subsection{Translational self-diffusion and rotational diffusion in pure water systems}
\brahms implements the diffusion measurement methods from Rapaport~\cite{rapa}~(pp.~124-128, 226).
\paragraph{{\tt diffusion}} ({\em logical}) Activates diffusion calculation.
\paragraph{{\tt nBuffDiffuse}} number of sets of data being collected at any time.
\paragraph{{\tt nValDiffuse}} Number of measurements contributing to the set used to produce a single unaveraged estimate of the diffusion coefficient.
\paragraph{{\tt stepDiffuse}} [number of steps] Interval between individual diffusion evaluations.
\paragraph{{\tt limitDiffuseAv}} Number of individual evaluations used to produce a single average which is then written on file.   
